##########################################################################################

library(data.table)
library(optparse)
library(ggplot2)
library(Seurat)
library(ggsci)
library(ArchR)
library(ggplotify)

##########################################################################################
option_list <- list(
    make_option(c("--rna_file"), type = "character"),
    make_option(c("--atac_file"), type = "character"),
    make_option(c("--scriptPath"), type = "character"),
    make_option(c("--out_path"), type = "character") 
)

if(1!=1){
    
    ## 单细胞表达文件
    rna_file <- "~/20231121_singleMuti/results/qc_atac/testis_combined.annotationCellType.qc.Rdata"

    ## atac文件
    atac_file <- "~/20231121_singleMuti/results/qc_atac/testis_combined_peak.combineRNA.qc.Rdata"

    ## atac原始的和rna取过交集的
    atac_raw_file <- "~/20231121_singleMuti/results/atac_res/testis_combined_peak.Rdata"

    ## 既往研究整理的代码
    scriptPath <- "~/20231121_singleMuti/scripts/scScalpChromatin"

    ## 输出
    out_path <- "~/20231121_singleMuti/results/celltype_plot/cell_qc"

}

###########################################################################################
parseobj <- OptionParser(option_list=option_list, usage = "usage: Rscript %prog [options]")
opt <- parse_args(parseobj)
print(opt)

rna_file <- opt$rna_file
atac_file <- opt$atac_file
scriptPath <- opt$scriptPath
out_path <- opt$out_path

dir.create( out_path , recursive = T)

###########################################################################################

a <- load(rna_file)
DefaultAssay(scrnat) <- "MAGIC_RNA"
DefaultAssay(scrnat) <- "RNA"
## scrnat

b <- load(atac_file)
## testis_combined_peak_combineRNA

##########################################################################################
## 已发表文献写好的脚本
source(paste0(scriptPath, "/plotting_config.R"))

###########################################################################################
## 样本类型颜色
col_sample <- c(
    rgb(red=212,green=31,blue=37,alpha=255,max=255) ,
    rgb(red=31,green=137,blue=66,alpha=255,max=255) ,
    rgb(red=39,green=45,blue=106,alpha=255,max=255) 
    )
names(col_sample) <- c("testis01" , "testis02" , "testis03")

width <- 4
height <- 3

###########################################################################################
# Violin plots of QC metrics

rna_ccd <- scrnat@meta.data

dodge_width <- 0.75
dodge <- position_dodge(width=dodge_width)

# scRNA nUMIs / cell
p <- (
    ggplot(rna_ccd, aes(x=batch, y=nCount_RNA, fill=batch))
    + geom_violin(aes(fill=batch), alpha=0.5, adjust = 1.0, scale='width', position=dodge)
    + geom_boxplot(alpha=0.5, width=0.25, outlier.shape = NA)
    + scale_color_manual(values=col_sample, limits=names(col_sample), name="batch", na.value="grey")
    + scale_fill_manual(values=col_sample)
    + guides(fill=guide_legend(title="batch"), 
      colour=guide_legend(override.aes = list(size=5)))
    + xlab("")
    + ylab("Number of UMIs per cell")
    + theme_BOR(border=FALSE)
    + theme(panel.grid.major=element_blank(), 
            panel.grid.minor= element_blank(), 
            plot.margin = unit(c(0.25,1,0.25,1), "cm"), 
            #aspect.ratio = aspectRatio, # What is the best aspect ratio for this chart?
            legend.position = "none", # Remove legend
            axis.text.x = element_text(angle = 90, hjust = 1))
    + scale_y_continuous(limits=c(0,50000), expand = c(0, 0))
)

out_file <- paste0(out_path, "/scRNA_nUMIs_per_cell_violin.pdf")
ggsave(  out_file , p ,width=width, height=height)

###########################################
# scRNA nGenes / cell
p <- (
    ggplot(rna_ccd, aes(x=batch, y=nFeature_RNA, fill=batch))
    + geom_violin(aes(fill=batch), alpha=0.5, adjust = 1.0, scale='width', position=dodge)
    + geom_boxplot(alpha=0.5, width=0.25, outlier.shape = NA)
    + scale_color_manual(values=col_sample, limits=names(col_sample), name="batch", na.value="grey")
    + scale_fill_manual(values=col_sample)
    + guides(fill=guide_legend(title="batch"), 
      colour=guide_legend(override.aes = list(size=5)))
    + xlab("")
    + ylab("Number of genes per cell")
    + theme_BOR(border=FALSE)
    + theme(panel.grid.major=element_blank(), 
            panel.grid.minor= element_blank(), 
            plot.margin = unit(c(0.25,1,0.25,1), "cm"), 
            #aspect.ratio = aspectRatio, # What is the best aspect ratio for this chart?
            legend.position = "none", # Remove legend
            axis.text.x = element_text(angle = 90, hjust = 1))
    + scale_y_continuous(limits=c(0,12000), expand = c(0, 0))
)

out_file <- paste0(out_path, "/scRNA_nGenes_per_cell_violin.pdf")
ggsave(  out_file , p ,width=width, height=height)

###########################################
# scRNA pctMito / cell
rna_ccd$percent.mt <- 100 * rna_ccd$percent.mito
p <- (
    ggplot(rna_ccd, aes(x=batch, y=percent.mt, fill=batch))
    + geom_violin(aes(fill=batch), alpha=0.5, adjust = 1.0, scale='width', position=dodge)
    + geom_boxplot(alpha=0.5, width=0.25, outlier.shape = NA)
    + scale_color_manual(values=col_sample, limits=names(col_sample), name="batch", na.value="grey")
    + scale_fill_manual(values=col_sample)
    + guides(fill=guide_legend(title="batch"), 
      colour=guide_legend(override.aes = list(size=5)))
    + xlab("")
    + ylab("Percent Mitochondrial Reads")
    + theme_BOR(border=FALSE)
    + theme(panel.grid.major=element_blank(), 
            panel.grid.minor= element_blank(), 
            plot.margin = unit(c(0.25,1,0.25,1), "cm"), 
            #aspect.ratio = aspectRatio, # What is the best aspect ratio for this chart?
            legend.position = "none", # Remove legend
            axis.text.x = element_text(angle = 90, hjust = 1))
    + scale_y_continuous(limits=c(0,50), expand = c(0, 0))
)

out_file <- paste0(out_path, "/scRNA_pctMito_per_cell_violin.pdf")
ggsave(  out_file , p ,width=width, height=height)


###########################################################################################

atac_ccd <- testis_combined_peak_combineRNA@cellColData %>% as.data.frame()

# scATAC TSS / cell
## 转录起始位点 (TSS) 富集分数- TSS 富集计算是信噪比计算。收集 TSS 参考集周围的读数，以形成以 TSS 为中心并在任一方向延伸至 1000 bp 的读数聚合分布（总共 2000 bp）。
## 然后，通过获取分布每个末端侧翼 100 bps 内的平均读取深度（总共 200 bp 的平均数据）并计算每个位置相对于该平均读取深度的倍数变化，对该分布进行归一化。
## 这意味着侧翼应从 1 开始，如果转录起始位点（基因组的高度开放区域）存在高读取信号，则信号应增加直至中间的峰值。
## 我们将归一化后分布中心的信号值作为我们的 TSS 富集度量。用于评估 ATAC-seq。
## TSSenrichment score： 也是文献中常见的一个细胞识别过滤参数，细胞保持活性需要一定数量基因的表达翻译，而基因转录时就需要开放 TSS 区域，
## 好结合转录因子、DNA 聚合酶等转录起始复合物，一般情况下 TSS 区域的 fragments 能占到细胞全部 fragments 的25%~35%，

p <- (
    ggplot(atac_ccd, aes(x=Sample, y=TSSEnrichment, fill=Sample))
    + geom_violin(aes(fill=Sample), alpha=0.5, adjust = 1.0, scale='width', position=dodge)
    + geom_boxplot(alpha=0.5, width=0.25, outlier.shape = NA)
    + scale_color_manual(values=col_sample, limits=names(col_sample), name="Sample", na.value="grey")
    + scale_fill_manual(values=col_sample)
    + guides(fill=guide_legend(title="Sample"), 
      colour=guide_legend(override.aes = list(size=5)))
    + xlab("")
    + ylab("TSS Enrichment")
    + theme_BOR(border=FALSE)
    + theme(panel.grid.major=element_blank(), 
            panel.grid.minor= element_blank(), 
            plot.margin = unit(c(0.25,1,0.25,1), "cm"), 
            #aspect.ratio = aspectRatio, # What is the best aspect ratio for this chart?
            legend.position = "none", # Remove legend
            axis.text.x = element_text(angle = 90, hjust = 1))
    + scale_y_continuous(limits=c(0,40), expand = c(0, 0))
)

out_file <- paste0(out_path, "/scATAC_TSS_per_cell_violin.pdf")
ggsave(  out_file , p ,width=width, height=height)


###########################################
# scATAC log10 nFrags / cell
atac_ccd$log10nFrags <- log10(atac_ccd$nFrags)
p <- (
    ggplot(atac_ccd, aes(x=Sample, y=log10nFrags, fill=Sample))
    + geom_violin(aes(fill=Sample), alpha=0.5, adjust = 1.0, scale='width', position=dodge)
    + geom_boxplot(alpha=0.5, width=0.25, outlier.shape = NA)
    + scale_color_manual(values=col_sample, limits=names(col_sample), name="Sample", na.value="grey")
    + scale_fill_manual(values=col_sample)
    + guides(fill=guide_legend(title="Sample"), 
      colour=guide_legend(override.aes = list(size=5)))
    + xlab("")
    + ylab("log10 nFrags")
    + theme_BOR(border=FALSE)
    + theme(panel.grid.major=element_blank(), 
            panel.grid.minor= element_blank(), 
            plot.margin = unit(c(0.25,1,0.25,1), "cm"), 
            #aspect.ratio = aspectRatio, # What is the best aspect ratio for this chart?
            legend.position = "none", # Remove legend
            axis.text.x = element_text(angle = 90, hjust = 1))
    + scale_y_continuous(limits=c(0,5), expand = c(0, 0))
)

out_file <- paste0(out_path, "/scATAC_log10nFrags_per_cell_violin.pdf")
ggsave(  out_file , p ,width=width, height=height)


###########################################
# scATAC FRIP / cell
## FRiP表示的是位于peak区域的reads的比例
p <- (
    ggplot(atac_ccd, aes(x=Sample, y=FRIP, fill=Sample))
    + geom_violin(aes(fill=Sample), alpha=0.5, adjust = 1.0, scale='width', position=dodge)
    + geom_boxplot(alpha=0.5, width=0.25, outlier.shape = NA)
    + scale_color_manual(values=col_sample, limits=names(col_sample), name="Sample", na.value="grey")
    + scale_fill_manual(values=col_sample)
    + guides(fill=guide_legend(title="Sample"), 
      colour=guide_legend(override.aes = list(size=5)))
    + xlab("")
    + ylab("log10 nFrags")
    + theme_BOR(border=FALSE)
    + theme(panel.grid.major=element_blank(), 
            panel.grid.minor= element_blank(), 
            plot.margin = unit(c(0.25,1,0.25,1), "cm"), 
            #aspect.ratio = aspectRatio, # What is the best aspect ratio for this chart?
            legend.position = "none", # Remove legend
            axis.text.x = element_text(angle = 90, hjust = 1))
    + scale_y_continuous(limits=c(0,1.0), expand = c(0, 0))
)

out_file <- paste0(out_path, "/scATAC_FRIP_per_cell_violin.pdf")
ggsave(  out_file , p ,width=width, height=height)


###########################################
###密度图
sampleNames <- unique(atac_ccd$Sample)
filterSample_cellNum <- table(atac_ccd$Sample)#过滤后样本数
sampleplot_list <- list()
proj.filter <- testis_combined_peak_combineRNA

for (i in 1:length(sampleNames)) {
  proj.i <- proj.filter[proj.filter$Sample == sampleNames[i]]#提取每个样本ArchRproject做密度图
  
  p <- ggPoint(
    x = log10(proj.i$nFrags),
    y = proj.i$TSSEnrichment,
    colorDensity = TRUE,
    continuousSet = "blueYellow",
    xlabel = "Log10(Unique Fragments)",
    ylabel = "TSS Enrichment")+ 
    geom_hline(yintercept = 4, lty = "dashed") + 
    geom_vline(xintercept = 3, lty = "dashed")+
    ggtitle(paste0(sampleNames[i],"\n","Cells Pass Filter =",filterSample_cellNum[[i]]))
  
  sampleplot_list[[i]] <- p
}

out_file <- paste0(out_path, "/scATAC_TSSenrich-Frag.pdf")

pdf(out_file)
print(sampleplot_list[[1]])
print(sampleplot_list[[2]])
print(sampleplot_list[[3]])
dev.off()

## 聚类图
p <- plotEmbedding(ArchRProj = testis_combined_peak_combineRNA, 
    colorBy = "cellColData", name = "TSSEnrichment", embedding = "UMAP", plotAs="points" , pal = paletteContinuous(set = "blueYellow"))
p <- as.ggplot(p)

out_file <- paste0(out_path, "/scATAC_TSSenrich.umap.pdf")
ggsave(  out_file , p ,width=width*2, height=height*3)

###########################################
## Fragment size distributions 
##  在ATAC-seq实验中，Fragment size distributions（片段大小分布）是指测量得到的DNA片段在长度上的分布情况。这对于理解DNA片段在基因组上的分布和ATAC-seq实验的质量至关重要。
## 通常，ATAC-seq的DNA片段大小分布在某种程度上应该符合预期的模式，但这也可能因实验设计、细胞类型和实验条件而异。以下是一些参考的期望模式：
## 1、Nucleosome-Free Region (NFR): 期望在片段大小分布中观察到清晰的NFR，即缺乏核小体的区域。这通常表现为片段大小在100-200 bp的峰。
##      缺乏核小体的区域通常代表着染色质的开放结构，这意味着在这些区域上DNA相对较为容易被核酸酶或其他切割酶切割，形成相对较短的DNA片段。
##      这种开放的染色质结构通常与基因的转录活动相关。
## 2、Mono-Nucleosome Peaks: 在片段大小为150-250 bp的范围内，你可能会观察到核小体区域的峰，这是由于单个核小体引起的。
## 3、Di-Nucleosome Peaks: 片段大小更大的峰可能表示双核小体区域，表明核小体紧密堆叠在一起。
## 4、Large Fragments: 大于500 bp的片段可能是由于染色质的结合或其他非特异性DNA断裂引起的。这可能是一些背景噪声，但也可能包含一些重要的生物学信息。
p <- plotFragmentSizes(ArchRProj = testis_combined_peak_combineRNA)
p <- as.ggplot(p)
out_file <- paste0(out_path, "/scATAC_FragmentSizes.pdf")
ggsave(  out_file , p ,width=width/2, height=height)

## TSS enrichment profiles
## ATAC-seq的TSS（转录起始位点）富集分布在解释基因表达和染色质结构方面具有重要意义。TSS富集分布的合理性可以通过比较样品之间的差异以及与已知生物学特征的一致性来评估。
## 一般而言，TSS富集分布的合理特征包括：
## 1、清晰的TSS峰值: 在TSS周围（通常是基因启动区域）应该存在清晰的ATAC-seq峰值，这反映了基因的转录活动。这些峰值通常表现为ATAC-seq片段的增加，提示染色质在这些区域上更加开放。
## 2、TSS附近的峰值对称性: ATAC-seq峰值在TSS附近通常是对称的，即在TSS两侧对称地存在ATAC-seq片段富集。
## 3、TSS富集与基因表达的一致性: TSS富集分布应与已知的基因表达模式一致。对于激活的基因，你应该观察到在其TSS附近的ATAC-seq峰值。
## 这种一致性可以通过比较ATAC-seq数据和RNA-seq等基因表达数据来验证。
## 4、启动子区域的清晰边界: 富集分布应该具有清晰的启动子区域边界，表明在这些区域上存在开放的染色质。
p <- plotTSSEnrichment(ArchRProj = testis_combined_peak_combineRNA)
p <- as.ggplot(p)
out_file <- paste0(out_path, "/scATAC_TSSEnrichment.pdf")
ggsave(  out_file , p ,width=width/2, height=height)

###########################################
## 每个细胞类型画基因以及counts的分布
width <- 10
height <- 9
#调整细胞顺序
grp_order = c("SSC", "Differenting&Differented SPG", "Leptotene",
"Zygotene", "Patchytene", "Diplotene",
"Early stage of spermatids", "Round&ElongateS.tids", "Sperm",
"Leydig cells", "Myoid cells", "Pericytes",
"Sertoli cells", "Endothelial cells", "NKT cells", "Macrophages")

rna_ccd <- scrnat@meta.data
rna_ccd$cell_type <- factor(rna_ccd$cell_type, levels = grp_order)

dodge_width <- 0.75
dodge <- position_dodge(width=dodge_width)

#细胞上色
use_colors <- c(pal_npg("nrc")(10) , pal_jco("default")(6))
#names(use_colors) <- unique(scrnat$cell_type)
names(use_colors) <-c( "Myoid cells","Leydig cells" ,          
"Endothelial cells","Zygotene",         
"Round&ElongateS.tids","Patchytene",
"SSC","Sperm" ,                 
"Diplotene","Early stage of spermatids", 
"Leptotene","Sertoli cells",            
"Macrophages","Differenting&Differented SPG",
"Pericytes","NKT cells")
##genecount/celltype统计基因数目并除以1000展示
rna_ccd$nFeature_RNAs <- rna_ccd$nFeature_RNA/1000
p1<- (
    ggplot(rna_ccd, aes(x=cell_type, y=nFeature_RNAs, fill=cell_type))
    + geom_violin(aes(fill=cell_type), alpha=0.5, adjust = 1.0, scale='width', position=dodge)
    + geom_boxplot(alpha=0.5, width=0.25, outlier.shape = NA)
    + scale_color_manual(values=use_colors, limits=names(use_colors), name="cell_type", na.value="grey")
    + scale_fill_manual(values=use_colors)
    + guides(fill=guide_legend(title="cell_type"), 
      colour=guide_legend(override.aes = list(size=5)))
    + xlab("")
    + ylab("gene number(x10³)")
    + theme_BOR(border=FALSE)
    + theme(panel.grid.major=element_blank(), 
            panel.grid.minor= element_blank(), 
            plot.margin = unit(c(0.25,1,0.25,1), "cm"), 
            #aspect.ratio = aspectRatio, # What is the best aspect ratio for this chart?
            legend.position = "none", # Remove legend
            axis.text.x = element_text(angle = 90, hjust = 1))
    + scale_y_continuous(limits=c(0,12), expand = c(0, 0))
)

out_file <- paste0(out_path, "/scRNA_nGenes_per_cell_violin.divideByCell.pdf")
ggsave(  out_file , p1,width=width, height=height)

#UMIs数目/cell
rna_ccd$log10nCount_RNA <- log10(rna_ccd$nCount_RNA)
p2<- (
    ggplot(rna_ccd, aes(x=cell_type, y=log10nCount_RNA, fill=cell_type))
    + geom_violin(aes(fill=cell_type), alpha=0.5, adjust = 1.0, scale='width', position=dodge)
    + geom_boxplot(alpha=0.5, width=0.25, outlier.shape = NA)
    + scale_color_manual(values=use_colors, limits=names(use_colors), name="cell_type", na.value="grey")
    + scale_fill_manual(values=use_colors)
    + guides(fill=guide_legend(title="cell_type"), 
      colour=guide_legend(override.aes = list(size=5)))
    + xlab("")
    + ylab("Normalized UMI Number(log-count)")
    + theme_BOR(border=FALSE)
    + theme(panel.grid.major=element_blank(), 
            panel.grid.minor= element_blank(), 
            plot.margin = unit(c(0.25,1,0.25,1), "cm"), 
            #aspect.ratio = aspectRatio, # What is the best aspect ratio for this chart?
            legend.position = "none", # Remove legend
            axis.text.x = element_text(angle = 90, hjust = 1))
    + scale_y_continuous(limits=c(2.5,5), expand = c(0, 0))
)

out_file <- paste0(out_path, "/scRNA_nUMIs_per_cell_violin.divideByCell.pdf")
ggsave(  out_file , p2,width=width, height=height)
#导出含genenumber和UMInumber的tsv文件
subset_df <- subset(rna_ccd, select = c(cell_type, nCount_RNA, nFeature_RNA, batch))
write.table(subset_df, file = paste0(out_path, "/cell_gene_UMI_number.divideByCell.tsv"), sep = "\t", row.names = FALSE)
